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I. INTRODUCTION 



One of the original testing grounds for QCD has been the study of scahng violations 
in structure functions in nucleon deep- inelastic scattering (DIS). Over the last 25 years a 
wealth of information has been accumulated on the unpolarized structure functions, Fi^2, 
covering some three orders of magnitude of the momentum transfer squared, Q^, and prob- 
ing the small Bjorken-x region down to x = Q"^ /2Mv ~ 10"'*, where v is the energy transfer 
to the nucleon at rest, and M the nucleon mass. The data at large {Q^ ^ 5 GeV^) 
show no deviation from the perturbative QCD (pQCD) predictions obtained from solutions 
of the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations More 
recently, accurate data have also been accumulating on the spin- dependent structure 
functions, gi^2, the dependence of which will be necessary to confront the pQCD expec- 
tations. 

The need for an accurate description of the dependence of both unpolarized and po- 
larized nucleon structure functions is severalfold. Firstly, in modeling hard processes other 
than DIS, the ingredients needed to calculate cross sections and event rates are parton dis- 
tribution functions at a certain scale Q^. Because the energies involved in purely hadronic 
reactions are often considerably larger than in DIS experiments, it is imperative that evolu- 
tion of the distributions (which are extracted primarily from DIS data) to the higher be 
as reliable as possible. 

At the other end of the energy spectrum, describing the behavior of structure functions 
at low is no less important. While the dependence of parton distributions is accessible 
to pQCD, their x dependence is of course not directly calculable. Often one tries to connect 
the non-perturbative content of parton distributions to valence quark models of QCD P|, 
which are presumed to be good approximations at ^ IGeV^ (for a recent summary 
see Ref. and references therein). To compare with data one usually evolves the leading 
twist component to the higher of the DIS experiments. However, because this procedure 
involves evolution through a region where the strong running coupling constant, a^, is no 
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longer small, one should at least take next-to-leading order (NLO) corrections into account 
before seriously testing the reliability of these calculations. 

In the spin-dependent sector, while the kinematics in experiments involving polarization 
are more limited than in unpolarized processes (average values in recent g\ measurements 
only vary between 2 and 10 GeV^), knowledge of the dependence of polarized structure 
functions is essential in connection with questions concerning the spin content of the nucleon. 
Aside from accurately correcting for the scale dependence of first moments of polarized 
distributions 0, such as in the Bjorken sum rule, one also needs to understand the effects of 

evolution upon the shape of the proton, and in particular neutron, gi structure functions 
themselves Q. While the leading order pQCD corrections to gi have been known for some 
time, the NLO effects have until now only been calculated for non-singlet (NS) distributions. 
Results recently published by Mertig and van Neerven for the anomalous dimensions of 
singlet operators at NLO have finally made it possible to compute the NLO corrections to 
the polarized singlet quark and gluon distributions as well. 

In this paper we analyze the effects of NLO corrections on polarized and unpolarized 
parton distributions, in both the flavor singlet and non-singlet channels, as a function of 
X and Q^. In particular, we concentrate on the evolution of the nucleon gi structure 
function, and its effect on the x dependence at small and intermediate x. To study the 
effects numerically, we develop an accurate evolution algorithm based on an inverse 
Mellin transform method. 

Often one encounters the need to evolve distributions whose moments are not able to 
be expressed in simple analytic forms such as those used in standard data fitting in Refs. 
|T0Hl2|. For example, most studies of structure functions in low-energy models of the nu- 
cleon, including recent attempts [|13| to describe the nucleon's non-perturbative sea within 
the pion-cloud model |T^, can only be expressed numerically, with the output needing to 



be parametrized before the standard evolution and fitting techniques can be applied. To 
facilitate evolution of both parametrized and calculated distributions, our method therefore 
provides for the possibility of evolving input distributions given only at a limited number of 



points in Bjorken x. 

The paper is organized as follows. In Section II we review the essential features of the 
DGLAP evolution equations at 0{a1), summarizing for convenience the relevant formulas for 
evolution of parton distributions and structure functions. More comprehensive discussions 
can be found for example in Refs. |[T5|-|T^]. A detailed description of the numerical method 



employed to perform the evolution is given in Section III. The numerical results, for both 
unpolarized and polarized DIS, are presented in Section IV, and conclusions are drawn in 
Section V. Finally, in the Appendices we collect all the relevant formulas for the LO and 
NLO anomalous dimensions, together with details about their extrapolation into complex 
moment space. 

II. EVOLUTION EQUATIONS AT NEXT-TO-LEADING ORDER 

In this Section we review the main features pertinent to the DGLAP evolution equations, 
including 0{a1) corrections. We begin with the usual forward Compton amplitude, T^^,, for 
the scattering of a virtual photon from a nucleon, 

T,u = ^jd'ze^^-{p,S\T{J,{z)JM)\P.S). (1) 

where p and q are the nucleon and virtual photon four-momenta, and S is the nucleon spin 
vector. The fact that in the Bjorken limit (p ■ g, —q^ oo) the amplitude T^i, is light-cone 
dominated allows for an operator product expansion (OPE) analysis of the time-ordered 
product T {J fj_{z) J , through which the long and short distance physics are formally 
separated. The singular bilocal operator T {J^{z)Jy{Q)) is expanded in a basis of regular 
local operators, with the light-cone singularity contained in the Wilson expansion coefficient 
functions [|l8H20l . 



The forward Compton amplitude T^y is related to the nucleon DIS hadronic tensor W^y 
through the optical theorem: 

W^y = -^mT^y. (2) 
vr 



The hadronic tensor is usually written in terms of the spin-independent Fi and F2, and 
spin-dependent gi and g2 structure functions as: 



where x = Q'^/2p-q and = —q^. The structure functions are related to parton distribution 
functions (light-cone momentum distributions) by convoluting the latter with the Wilson co- 



efficient functions [0,^,^], describing the hard (perturbatively calculable) photon-parton 



interaction. For the proton's F2 structure function, for instance, one has: 



F^'ix, Q') = f [C^' [x, Q'^} xq^' Q ^ 

+ lci{x,Q') ^^[^,Q')+ lc^{^,Q') ^^(^'^'))' (4) 

where q^^{x, Q'^) = {u + u — d — d — s — s){x, Q"^) is the flavor non-singlet combination (for 
three flavors), while the singlet component is S(x, Q^) = J2giQ + q){^^Q'^)^ and G{x,Q'^) is 
the gluon distribution. 

Similarly for polarized DIS, one has for the gi structure function of the proton: 



Jx V \ 



ix y 

+ (x, Q') AS (^^, Q'^ + ^-AC^ (x, Q') AG (^^, ) , (5) 

where the polarized singlet distribution is AS(x, Q^) = I]g(Ag + Aq){x,Q^), AG{x,Q^) is 
the polarized gluon distribution, and the non-singlet combinations Aq^ and Aqg are (for 
three flavors) given by: 

Ag3(x, Q2) = (Am + Au- Ad- Arf) (x, Q^), (6a) 
Ag8(x, Q2) = (Am + Am + Ad + Arf - 2 (As + As)) (x, Q^). (6b) 

Note that the polarized quark singlet coefficient function AC'^{x,Q'^) = 5(1 — x) + 0(«s), 
while the polarized gluon enters only at NLO, AC"^(x, Q"^) = 0{as)- We shall return to this 
point in Section |IVB| . 



A. Evolution of Parton Distributions 



The evolution of parton distribution functions is most naturally expressed in terms 
of their moments, by solving the renormalization group equations. Defining the nth Mellin 
moment of a parton distribution g(x, Q^) by: 
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one obtains for the NS moments: 

^2\ ^ f ,,2\ n 



(7) 
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(8) 



^ 2/3o 2(31 

where 7/5*5'''" are the anomalous dimensions at leading and next-to- leading orders, and /x^ is 
some reference scale. The constants /3o, I3i are the expansion coefficients of the Gell-Mann- 
Low function I3{Q'^) up to next-to-leading order: 



(3o = 11~\nj, = 102 - yiV;, 



(9) 



for Nf active quark flavors. Note that Eq.(|^) is expressed within the MS factorization scheme 
— for the connection to other factorization schemes see Appendix |B[ 

Because of mixing between local operators of the OPE which are singlets under SU(iV/) 
transformations, one has coupled evolution equations for the singlet parton densities: 
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where R", P!^ and A± are given by: 
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For polarized parton distributions one has identical evolution equations to those 
in Eqs.(H) and (p!OD, with Qn,S„,^„ replaced by their spin-dependent counterparts 
AQ„,AS„,A^„. The differences in the evolution arise from the anomalous dimensions, 
which are in general different for polarized and unpolarized processes. The explicit expres- 
sions for these at LO and NLO are given in Appendix 0. 

The anomalous dimensions are also related to the Mellin moments of the splitting func- 
tions Pij{z) = q,G): 
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(12) 



where Pij{z) gives the light-cone probability of finding a parton i inside parton j, carrying 
a fraction z of the parent parton's light-cone momentum. If the splitting functions (or the 
anomalous dimensions) are known, then once the input parton distribution is fixed at a 
scale /i^, the dependence of its moments can be calculated via Eqs.(||) and (p^Ol). The x- 
dependent distribution at the new scale is then obtained by inverting the Mellin transform 
in Eq.(0). 



B. Evolution of Structure Functions 

From the evolved parton distributions one can obtain the measured structure functions 
in two ways. Firstly, one can evolve the parton distributions, then convolute the result 
with the corresponding Wilson coefficient functions, as in Eqs.(|[) and (|]). In practice, 
however, to accurately determine the convolution integral requires interpolating the evolved 
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distribution, thereby introducing an additional source of numerical error. Alternatively, the 
evolved moments can be multiplied by the corresponding moments of the Wilson coefficient 
functions, before reconstructing the structure function via the inverse Mellin transform. In 
this work we adopt this second method to evaluate the structure function evolution. 

The explicit relations between the moments of the non-singlet and singlet parts of the 
F2 structure function and the corresponding moments of the parton distributions are: 



(13) 
(14) 



For the F3 structure function, in the case of neutrino-nucleon scattering, one has: 



(15) 



where Q^{Q'^) represents the nth moment of the total valence distribution. Up to NLO in 
the MS scheme the Wilson coefficients are given by: 
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with Cp = 4/3 and Tp = Nf/2. The functions 5*1 and 5*2 are given in Appendix [A 1| . 

The analogous expression for the moments of the polarized structure function gi{x, Q'^) 
(for Nf = 3) is: 



(19) 



where ± refers to the proton and neutron, respectively. The polarized Wilson coefficients 



ACn are given by P,23 



Acr (Q^) = ^c^XQ') = AC'f + ^4?^AC(^)^ 
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As in the unpolarized case, one has to use anomalous dimensions and Wilson coefficients 
calculated consistently within the same factorization scheme. For the NLO evolution of 
both unpolarized and polarized distributions we work within the MS-scheme throughout. A 
property of this scheme is that AC^^i = 0, resulting in zero gluonic contribution to the first 
moment of gi, and hence no correction to the Ellis- Jaffe sum rule [03]. All other moments 
of the structure function gi{x^Q'^) do, however, contain gluonic corrections. 



III. NUMERICAL SOLUTION 



Having outlined the essential features of the DGLAP evolution equations, we next ex- 
amine methods which can be used to efficiently and reliably extract their solutions. The 
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most direct approach involves simply integrating the integro-differential version of the equa- 
tions directly (see e.g. Ref. [^]). With this method, to obtain evolved distributions in the 
range [xq, 1] the input distribution needs to be known only for x > Xq- However, many 
iterations are necessary to obtain accurate results when evolving over a large range of Q^, 
which consequently decreases both the efficiency and accuracy of the evolution. Other pro- 
posed methods have made use of complete sets of orthogonal polynomials, such as Bernstein 
polynomials ||2^ or Laguerre polynomials |2^. Difficulties exist, however, in applying these 
methods to regions of very small or very large x. Instead, we will follow an approach using 
the Mellin transform technique similar to that discussed in Ref. [0, but differing in a 
number of aspects which we now discuss. 



A. Contour Method 

The advantage of the Mellin transform (or contour integration) technique lies in its 
simplicity, and efficiency with computing time. In the past this method has been applied in 
cases where the moments of the input distributions were known in analytic form, such as 
when q{x) is a linear combination of terms of the form ax''{l — xY or similar [|12|. The need 
often exists, however, to evolve input distributions which are not able to be expressed in this 
form. For example, one may wish to evolve distributions obtained from model calculations 
0, which generally do not have a simple analytic form. Therefore it is necessary that 
one has available an accurate NLO evolution procedure for applications where the input 
distributions are known only at a limited number of points in x. 

The main steps involved in the Mellin transform technique can be summarized as follows: 

• The Mellin moments of the input distributions Qn(/U^) are calculated numerically at 
the input scale /i^. 

• The calculated moments are evolved to the scale via Eqs.(|^)and ([TI1|). 

• The parton distributions are reconstructed via the inverse Mellin transformation: 
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1 PC+lOO 

q{x,Q^) = TT- / dn x-'^QniQ^). (23a) 

ZTTl Jc-ioo 

In practice it is more convenient to express Eg. ( P3a| ) as an inverse Laplace transfor- 
mation: 

1 fC+ioo 
ZTTt Jc—ioo 

One disadvantage in principle of this method is that in order to calculate the moments 
of distributions, the input must be known in the entire x-region [0, 1]. However, as will be 
seen in Section the restriction to the interval [xq, 1] has negligible effects numerically. 

B. Evaluation of Moments 

If the input distribution is of the form g(x,/i^) = ax^{l — xY, its moments are simply 
expressed in terms of Beta functions B: 

If the analytic moments are not available, each Mellin moment must be integrated numer- 
ically. This causes difficulties, however, for complex values of n, since the integrand in 
Eg. ( |23 bp oscillates rapidly as x — > 0. This can be seen by writing n = c + id, where c, d are 
real: 

Q^(^^) = [' dx x-i+V^i"" qix, fj,^). (25a) 
Jo 

In practice it is easier to perform the integration by making the substitution Inx —>■ t, which 
gives an integrand that oscillates with constant freguency: 

Q„(^2)=r e(^-^)V''*g(e*,/i2). (25b) 

J —oo 

The real and imaginary parts of this integral are then evaluated separately by using a suitable 
adaptive guadrature algorithm (such as that provided by the IMSL numerical library routine 
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DQDAWO p9[). This way accurate moments can be obtained for any input distribution 
given as a function of Bjorken-x over the entire region, < x < 1. 

Two additional problems arise when using input that is given only at a limited number of 
x-points. Firstly, without an adaptive quadrature routine the moments are often inaccurate, 
with the error increasing with the imaginary part of n. One is then forced to interpolate 
the input between the x-points, which leads to a second problem, namely that the x-range 
covered by these input points is limited to xq < x < xi. 

Clearly the quality of the moments, and hence the accuracy of the evolved parton dis- 
tributions, depends strongly on the quality of the interpolation routine. For this we have 
used cubic splines and higher order B-splines, both with a 'not-a-knot' condition (see e.g. 
IMSL routines DCSINT and DBSINT, respectively p9[). Taking equidistant x-points leads 
to very good results everywhere except in the neighborhood of the boundaries at x = Xq and 
Xi. To obtain sufficiently good results over the entire x-range one should choose the location 
of x-points suitably, in particular more points near the boundaries than in between. 



C. Numerical Inversion 

The final step, which is numerically the most difficult, is the contour integration in the 
inverse Mellin (or Laplace) transformation. Depending on the value of Bjorken-x at which 
the parton distribution is to be reconstructed, two different methods are employed for the 
regions x < 0.80 and x > 0.80. 



In the X < 0.80 region Crump's method |Q is applied (using e.g. the IMSL routine 



DINLAP [^), in which the inverse Laplace transform of the function F{n) is approximated 
by a partial sum of a Fourier representation: 



fit) ^ - 

T 



1 °° f kilt kut 
-F(c) + Y me{F{c + kiri/r)) cos ^m{F{c + kiri/r)) sin 



(26) 



The value of c has to be on the right side of all singularities appearing in F{n). From 
Regge physics one obtains the conditions 3?e n > 2 and 3?e n > 1 for the singlet and 
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non-singlet distributions, respectively. The parameter r controls the step length in the 
imaginary direction. This method fails, however, when f{t) reaches values close to zero, 
which occurs for large x. In the kinematic range x > 0.8 we therefore utilize an alternative 
algorithm, in which the inverse Laplace transform is computed via a Laguerre expansion 
(see e.g. IMSL routine DSINLP |^|). The performance of this procedure has been checked 
by first performing the Mellin transform by computing the moments, and then transforming 
back with the above routines to compare the output with the original input. 

One further comment should be made about the small-x region. It is in general more 
problematic to compute the inverse Laplace transform f{t) numerically when t is very large, 
or equivalently when x is very small. In practice we find that the above procedure gives reli- 
able results for parton distributions at least down to x values reached by current experiments 
at HERA, namely x ~ 10"'^. 



D. Computation of the Running Coupling Constant 

If one performs evolution at low values of Q^, care must be taken when calculating 
the running coupling constant a^. At NLO, can be computed exactly by solving the 
transcendental equation numerically: 



KcD Po(^s Po 



47r ^ A 



0, (27) 



^QCD 

as obtained from the renormalization group analysis. For small values of as one can find 
approximate solutions by using the formula: 

"^^^^ Polr^QV^lcD PI Ih^QVA^cd' ^ ^ 

In practice this approximation turns out to be accurate only for ^ 1 GeV^. At smaller 
scales, ^ 0.4 GeV^, the error in as is about 7% (with Aqcd = 200 MeV in both Eq.(p7D 
and (^) and Nf = 3). In the region ^ 0.2 GeV^ from which low energy model calculated 
input is often evolved, use of the approximation Eq.(P^) can give errors of up to 20% in the 
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running coupling for the same values of Aqcd and Nf. The discrepancy becomes bigger 
still for larger Aqcd- Of course Aqcd itself can be extracted using Eq.(^) from data at 
large Q^, since there the difference between Eqs.(^) and ( pSj ) is minimal, however the use 
of Eq.(p8D outside this region could be problematic given the above differences. 



IV. RESULTS 

Using the numerical method described in Section III, we present here the results of 
the evolution of parton densities including all the NLO corrections. We should stress that 
our aim is not to perform a complete fit to the available data — such programs of data 



parametrizations are adequately addressed in Refs. [p!0|-[T^, for example. Our purpose will 
be to simply demonstrate the potential differences in the evolved distributions when order 
corrections are taken into account. After first establishing the degree of accuracy and 
the region of reliability of the evolution program in the unpolarized sector, we then turn our 
attention to the more topical discussion of evolution of polarized distributions at NLO. 



A. Unpolarized Distributions 

To demonstrate the effect of the numerical algorithm on a typical parton distribution, for 
the spin-averaged distributions we take as input the latest parametrization from the CTEQ 



collaboration 1 1 



xuvix) = 1.37x°-"^^(l - xf-'\l + 6.25x°-««°), (29a) 

xdvix) = 0.801x°-^^^(l - x^ '^l + 1.69x°-3^^), (29b) 

xG{x) = 0.738x-°-2^^(l - x)^-^\l + 7.30x), (29c) 

x{d{x) + u{x)) = 0.1094a;-°-2*^'^(l - xf-^\l + 17.5x), (29d) 

given at = 2.56 GeV^. Evolution at NLO is carried out in the MS scheme, using 
Aqcd = 239 MeV for Nf = 4 flavors, and taking into account the charm threshold at 
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=(1.6 GeV)^. (Although the effect of neglecting the charm threshold increases with 
decreasing Q^, even when evolving down to ~ 0.5 GeV^ it is still very small.) 

To illustrate the performance of the evolution algorithm, we first examine NLO evolution 
of the total valence distribution, x{uy + dy)- In Fig.l the program's accuracy is illustrated 
when using different numbers of points, N^., in Bjorken-x, with the input evolved to = 100 
GeV^. With = 50 points in the range [0.005,1], the results are indistinguishable from 
those obtained by calculating the moments of (|29a| , |29b|) analytically. In this case the effect of 
changing the lower limit of integration from xq = to Xq > is negligible. Excellent results 
can be attained for nearly all x even with A^^; = 25 (solid curve, marked "^25" in Fig.l). 
Decreasing the number of points to = 13 has very little effect in the intermediate-x region, 
0.2 ^ X ^ 0.8, but at larger x the accuracy becomes somewhat worse. For comparison, we 
also show in Fig.l (dashed curves) the results obtained by calculating moments via a simple 
histogram method, 

Qn{Q^)= I'dx x"-2 {xq{x,Q')) 

« 5: x,g(x„ Q') p ^_ ^~' ) ■ (30) 

For Nx = 50 the discrepancy with the interpolation method is rather more dramatic, partic- 
ularly at large x, and to obtain a similar degree of accuracy one needs A^^ in excess of 1000. 
Similar accuracy will also be obtained when evolving different input distributions, such as 
those from Refs. p,|10],|12 . 



In Fig. 2 we compare the valence x{uv+dv) distribution evolved with leading and next-to- 
leading order corrections. For the LO evolution we use Aqcd = 177 MeV, according to the 
latest fit results from Ref. ||1T|. To avoid ambiguities associated with the factorization scheme 
dependence of the (generally unphysical) NLO distributions, we work in the DIS scheme. 



(Differences between distributions calculated in different schemes are of order as{Q ) ^1 
and while negligible at high Q^, they can become quite large at ^ 1 GeV^.) Although 
in general input distributions at LO and NLO need not be identical, especially when one 
seeks precisions fits to DIS data, we will simply demonstrate the effects by assuming the 
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input shapes to be the same. This assumption is an imphcit one, for example, in low-energy 
effective model calculations of leading twist parton distributions where because the 

input model scale is a priori unknown, the same calculated distribution must be evolved 
with either LO or NLO corrections. 

The main effect of NLO corrections on the valence distributions (with upwards evolution) 
is to make them softer compared to those evolved with LO corrections only. Conversely, 
NLO evolution produces somewhat harder distributions when evolved downwards to low 
scales = 0.5 GeV^, such as those associated with valence distributions calculated from 
low-energy (constituent quark) models of QCD. 

In connection with this, we should note that uncertainty in the precise value of Aqcd 



(the Particle Data Group gives 234 ±56 MeV [^) can be translated into quite a significant 
difference in the starting scale Ql for the evolution of the calculated model results. Taking 
the upper and lower values of Aqcd = 290 and 178 MeV, one obtains essentially identical 
shapes of the evolved distributions with Qq = 1 and 0.5 GeV^, respectively. Even given 
the fact that the exact scale for such model calculations is a priori unknown, a factor 2 (!) 
represents quite a large effect in the phenomenology of the scale dependence of low-energy 
model distributions. 

Turning now to the singlet sector, we illustrate in Fig. 3 the effects of NLO corrections on 
the distribution T,{x,Q^). Taking the CTEQ parametrization at Ql = 2.56 GeV^ as input, 
the distribution is evolved to = 10 and 100 GeV^ using the DIS scheme to compare 
physical quantities at LO and NLO (in the DIS scheme xS(x,Q^) represents the structure 
function F2 measured in deep- inelastic neutrino- nucleon scattering). The results here have 
been obtained with A^^. = 60 on a logarithmic scale in Bjorken-x, and agree to within 0.01% 
with those using analytic moments. Similar accuracy is obtained also for the singlet gluon 
distribution. 

Having established the quality of performance of our numerical method for both non- 
singlet and singlet evolution of unpolarized distributions, in the next Section we discuss the 
NLO effects on polarized parton densities. 
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B. Polarized Distributions 



Evolution of polarized parton distributions in the non-singlet sector is straightforward, 
since at NLO it is governed by exactly the same anomalous dimensions as for unpolarized 
DIS. Therefore the effects are identical to those in Figs.l and 2. In the singlet sector, on 
the other hand, the anomalous dimensions for polarized and unpolarized scattering are very 
different. In fact, due to the complexity of the problem, the relevant splitting functions 
and anomalous dimensions necessary to describe the NLO singlet evolution as a function 



of X have only recently been calculated 0. A recalculation by Vogelsang |3^, using an 
independent method, has confirmed the results in |p. To use the inverse Mellin method for 
the evolution, it is necessary to analytically continue the results of Ref. |^ into the complex- 
n plane. In doing so care must be taken since, unlike for the unpolarized distributions, the 
continuations for spin-dependent structure functions are performed from odd moments. The 
full details of the continuation can be found in Appendix |A 2| 

To illustrate the role of NLO corrections in spin-dependent distributions, we have used 
as input the global fit by Gehrmann and Stirling (Set A) p6|: 



xAuvix) = 0.327a;°-^^(l - a;)^-''^(l + 18.36a;), (31a) 

xAdvix) = -0.139a;°-^^(l - x^'^^l + 18.36x), (31b) 

xAu{x) = xAd{x) = xAs{x) = xAc{x) = 0, (31c) 

xAG{x) = 16.64 a; (1 - a;)^■^^ (31d) 



with a starting scale of = A GeV^. Following Ref. |36|, in LO we again perform the 



evolution with three active fiavors and Ag^^ = 177 MeV, as for the unpolarized evolution 



^While completing this paper, similar NLO analyses of polarized distributions in Refs. |33.34| 
came to our attention. A different prescription for continuation of the singlet polarized anomalous 



dimensions into complex-n space is used in Ref. [33|, although we have checked that numerically 



these are in fact equivalent — see Appendix A 2 
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above. However, contrary to the procedure adopted in we omit the anomalous gluon 



term when calculating the structure functions gi{x,Q'^) and g'^{x,Q'^) at leading order. 
This term arises only at NLO, and in a consistent treatment should only be included in 
NLO evolution — see Eqs.(^) and (0). In NLO we also use three active flavors, and a 
different value of ^qcd = 239 MeV, corresponding to that obtained by CTEQ in their flts 
to unpolarized parton distributions. As in the discussion of the unpolarized evolution in 
Sec. IV A, for illustration purposes we take the same input shape for both LO and NLO 
evolution. 

In Fig. 4 the x dependence of AS(a;,(5^), evolved down to 2 GeV^ and up to 10 GeV^, 
is shown between x = 10~^ and 1. The difference between distributions evolved in LO and 
NLO can be sizeable below x ~ 5 x 10~^, and even larger with decreasing x. Furthermore, 
it is interesting to observe that at = 10 GeV^ even the curvature changes sign compared 
to the LO result. Similarly large differences between LO and NLO evolution appear for the 
polarized gluon distribution, AG{x, Q^), at small x — see Fig.5. These dramatic effects are 
largely independent of the shapes of the input LO and NLO distributions, and can be traced 
back to the dominant role played by NLO corrections at small x. This becomes evident if one 
considers the x — >• behavior of the splitting functions Pij{x) = P^j\x) + {as/4:T[) pj;j\x) 
PJ37[| — in LO one has: 





-0)- 


-4Cp, 


(32a) 




-0)- 


^ -8Tp, 


(32b) 


p!.%- 


-0)- 


8Cf, 


(32c) 


Pcci^ " 


-0)- 


- IQCa, 


(32d) 



while at NLO: 



Pil\^ ^ 0) ^ (sCfCa - IGCfTf - 12Cj) In^ x, (32e) 



P^G (a: ^ 0) ^ - (IQCaTf + SCfTf) In^ x, (32f ) 

Pen (x ^ 0) ^ (iGCaCf + 8Cj) In^ x, (32g) 



Pckx ^ 0) ^ (32Cl - IQCfTf) In^ x. (32h) 

For small enough x, the NLO corrections Pij\x) can become quite large and overwhelm 
the a<(/47r factor. Note also that even the sign of P^g\x) changes compared with Pjj^\x) as 
x^O (for Cf = 4/3, Tf = Nf/2, Ca = S and Nf>2). 

With the evolution of the polarized quark singlet and gluon distributions known, one 
can now examine the effect of NLO corrections on the spin-dependent structure function 
gi{x, Q^). In Fig.6 we show the proton structure function xgf^x, Q^), evolved to the averaged 
scale Q"^ = 10 GeV^ at which the SMC data were taken. The important result to notice 
is that the NLO corrections are large in the range of x covered by the experiment. In the 
region x ~ 0.2 the corrections can in fact be of the same order of magnitude as the quoted 
error bars on the data. Furthermore, the fact that these are negative is consistent with the 
known reduction of the first moment of gi [0 when NLO corrections are included. 

In Fig. 7 we also show the polarized neutron structure function xg'^{x,Q'^), evolved at 
LO and NLO down to a scale = 2 GeV^ for comparison with the SLAG E142 data 
0]. Because evolution here involves a somewhat lower resolution scale, the effects of NLO 
corrections should be be relatively larger. Glearly, one sees that these can have a potentially 
significant effect on the shape of the structure function, when the data points are evolved 
to the same value of Q^. 

V. CONCLUSION 

We have examined in detail the role of NLO corrections on the shapes of the nucleon's 
singlet and non-singlet parton distributions, using a evolution algorithm based on the 
inverse Mellin transform in complex moment space. For a complete description of polarized 
and unpolarized evolution at NLO it has been necessary to extrapolate recently obtained 
results for the polarized singlet NLO anomalous dimensions into the complex-n plane. Of 
particular importance are the corrections to the spin-dependent gi structure functions 
of the proton and neutron. We find that in the region of x and relevant to current 
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and upcoming experiments, the NLO effects are not negligible, and should be included 
in future analyses of data. Finally, we should make a note about the behavior of the 
polarization asymmetry Ai = g\jF\^ which was assumed to be independent of in previous 
data analyses P,^. While the non-singlet parts of g\ and F\ evolve similarly at NLO, the 
evolution of the singlet components is dramatically different, as seen in Figs. 3-5, especially 
at low X. Therefore one should be wary of the possible errors introduced into g\ at small x 
if one extracts this from the A\ data under this assumption. 
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APPENDIX A: ANOMALOUS DIMENSIONS 



1. Unpolarized DIS 



Although the results for the full anomalous dimensions up to O(a^) "^^^ be found through- 
out the literature, for convenience we catalogue them here in full. In the MS scheme the 
anomalous dimensions are given by the perturbative expansion: 



(0),n , f CisiQ 



2\\ 2 



4:71 



Arc 

(0),n ( 

^ I 47r 



(1)." , 
Ins +■ 



7 



{k),n 



(k'],n (k),n 

yrcq tgg j 



(Ala) 
(Alb) 

(Ale) 



where the order anomalous dimensions are P,|38 



AO),n _ (0),n _ 
Ins — Iqq — ^ ^' 



4^1 (n) -3 



n{n + 1) 



IqG 



(0),n 
TCq 



-8Tf 
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(0),n _ 9^ 



+ n + 2 
n{n + l){n + 2)' 
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45i(n) - ^ 



4 



3 ' 



(A2a) 
(A2b) 
(A2c) 
(A2d) 



n{n-l) (r2 + l)(n + 2) 
with Cp = 4/3, Tp = Nf/2 and Ca = 3. The anomalous dimensions at 0{a1) are given by 



3§: 



7jv5 — '-^i; 



lQSi{n){2n + I) 
n^in + 1)2 



+ 16 2Si{n) 



n{n + 1)^ 



{S,{n)~S',{n/2)) 



+24:S2{n) + 645(n) - 85^(^/2) - 3 



3n3 + - 1 



+CaCf 



n^{n + 1)^ 
536 



16?7 



2n^ + 2n + 1 



9 



n^(n + 1)^ 
^i(n) - 8 ( 25i(n; 



— S2(n)-325(n)+45^(n/2)-- 



n{n + 1)^ 
17 



(252(n) - S^(n/2)) 
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4 151n^ + 236n^ + 8871^ + 3n + 18 + 2n + 1 



160^,, 32^,, 4 1611n2 + 5n-3 



+ 32n^ + 49n^ + 38n^ + 28n + 8 
(n-l)n3(n + l)3(n + 2)2 ' 



IqG 



—SCaTf 
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Tgg 



CaTf 
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+CfTf 



8 + 16 



(n - l)n3(n + + 2) 

2n^ + 5n^ + Sn^ + 7^2-2^-2 64 

~ y 



5*1 (n) + 645*1 (n) , ,^ ^, 

9 ^ ^ ^ ^ (n-l)2n2(^ + 1)2(^ + 2)2 



n +n + l 



+32^'(n/2) — 

^ ^(r2-l)n(n + l)(n + 2) 

4 457^9 + 2742n^ + 6040n^ + 6098n6 + 1567^^ - 2344ra^ - 1632^3 

~9 (n- l)2n3(n + l)3(r2 + 2)3 

560n2 + 1488n + 576 



(n- l)2n3(n + l)3(n + 2)2 



16Si{n)S'^{n/2) + 32S{n) - 4S^(n/2) 



(A3e) 



The expressions (^^) and ( [A3D are also valid for complex n if the following analytic 
continuations are used [IT^l: 



where 



n , 

j=iJ 

n 1 

S2{n) = E-^-C(2)-vi/'(n + l), 

^sH = E--C(3) + nn + i), 
j=ij 



1(1 + v)Si{n/2) + 1(1 - r^)^K(^ - l)/2), 1 = 1,2, 3, 

-ly 



i=i 



-k(3) + ^ 1^ - ^ [^{{n + l)/2) - M/(n/2)] 



Jo 1 + x J 



(A4a) 
(A4b) 
(A4c) 



(A4d) 



(A4e) 



7^ = 0.577216, 
C(2) = 7rV6, 
C(3) = 1.202057, 



(A5a) 
(A5b) 
(A5c) 



and the Polygamma functions, are defined as 
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The integral over the Dilogarithm \a2{x) in S{n) can be calculated analytically by using the 
approximation: 
Li2(x) 



1 + x 



0.0030 + 1.0990X - 1.5463x^ + 3.2860x'^ - 3.7887x^ + 1.7646x^ (A7) 



obtained from a least squares fit. Similar to the fit in Ref. [0, it is a polynomial of degree 
five, but gives slightly better values for the anomalous dimensions, as can be seen by checking 
the following relations for the n = 2 moments: 

1^^' = (A8a) 

if^'""^' = -1^^'- (A8b) 

The Polygamma functions must be computed numerically. For the case \n\ > 10 one 
applies the asymptotic expansions: 

. s , 1 1 1 1 /» X 

\E'(n)~lnn -, (A9a) 

^ ^ 2n 12^2 120^4 256^6' ^ ^ 

■^'(n) ~ - H r H H (A9b) 

^ ' n 2n2 6n3 30ri5 42^^ 30^9' ^ ' 

^^"(n] ^— + — + — — fA9cl 

^ '~ n2 n3 2n^ Qn^ 6^8 lO^io Qn^^' ^ ' 

To obtain greater accuracy for \l'"(n) we apply an Euler transformation to the asymptotic 

series in ( |A9c| ) , ending up with: 

„ 11 1 1 1 3 5 

^^''~~^~^~2^^6^~T6^~20^~48^" ^ ' 

In the case of |n| < 10 the recursion formula 

^M(ri + 1) = m^"'\n) + (-l)"^m! n^"™ (AlO) 



is utilized repeatedly |12| . 

At next-to-leading order, a factor (—1)" appears in the anomalous dimensions [^0 
(Eg. ( |A4| )). Therefore the analytic continuation depends on whether the evolution equa- 
tions for a particular parton distribution are valid for even or odd moments, in which case: 
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-1^ 



(All) 



for n even or odd, respectively. The 1] factors relevant for various combinations of parton 
distributions can be derived from the crossing relations: 



F^{-x,Q^) = +F^{x,Q^), 
i^2,3(-a:,Q') = -F2,3(x,Q2 



(Al2a) 
(A12b) 



For the non-singlet combinations these are (for Nf = 4 active flavors) ||T2| 

u — u ; ?7 = — 1 

d — d -,1] = —1 

{u + u) — {d + d) ; = +1 

{u + u) + {d + d)-2{s + s) -,7] = +1 

{u + u) + {d + d) + s + s-3{c + c) ;?7 = +l. 

The singlet distributions are evolved with rj = +1. 

2. Polarized DIS 



(A13a) 
(Al3b) 
(Al3c) 
(AlSd) 
(Al3e) 



For completeness, we also list the full anomalous dimensions relevant for polarized dis- 
tributions at NLO. The order anomalous dimensions are given by [@,^: 
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Concerning the order a^, the anomalous dimension 7^5 for the non-singlet operator is the 
same as in the unpolarized case. The singlet dimensions, as evaluated in Ref. [H, are as 



-Tf. 
3 



(A14a) 
(A14b) 
(Al4c) 
(Al4d) 
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follows 
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Note the correction to the T^g'^jTc^'" anomalous dimensions in Q. 
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n+1 ■ (n+l)2 (n + l)3 

To make use of the anomalous dimensions in the contour evolution method, they must 
be analytically continued into the complex-n plane. An important fact to observe is that 
due to the crossing relation: 



gi{-x,Q'^) = -gi{x,Q'^) 



(A16) 



the OPE relates only odd moments of the structure function gi^x^Q"^) to the sum of pos- 
sible twist-two operator matrix elements between nucleon states (each multiplied by the 
appropriate moments of the Wilson coefficient functions). This is contrary to the case of 
F2(x, (5^), where an analogous relation holds for n even. To evolve AE(a;, Q^), Ag3(x, Q^) 
and Ag8(x, Q"^) one therefore has to continue the anomalous dimensions from odd n. 

In addition to the functions 5'i 2,3, ^ unpolarized anomalous dimensions, 

for polarized DIS one must also determine the analytic continuation of the functions: 



j=i J 

j=l J 

j=l J 



(A17) 
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Specifically, this must be done for the case of 5"!, 5*2, 5*3, S'1^2, 'S'2,1, and 5'i_2- Because all 
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these functions appear with argument (n — 1) in Eq.( |A15|) (see below), one must continue 
them into complex n starting from even integers n. Using the relation: 
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with ?7 = 1) for n even (odd), a straightforward calculation leads to: 
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for n even. The equality 111 
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then gives: 
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The function S'i^2('^) can be obtained from the above expressions via |]39|j41[| : 



Si^2{n) = Si{n)S2{n) + S^^n) - S2,i{n). 



(A24) 



The most difficult function to continue is 5*1,2 (^)- For this we start with (see Eq. (|A2CI|)) 
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After some manipulation one finds: 
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(A27) 



where 2-^i('^, c, z) is the hypergeometric function. For our apphcation it is convenient to 
use the integral representation for 2F1 ||4^ : 

Tic) rl 



2Fi{a,b,c, z) 



dtf-\i-ty-'-\i-tz)-\ 



T{b)T{c - b) Jo 
Finally, making the substitution t ^ one arrives at: 

2Fi(l, n+l,n + 2,z)+ 2^1(1, n + l,n + 2,-z) = 2(n + 1) / " du 
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APPENDIX B: FACTORIZATION SCHEME DEPENDENCE 

In this Appendix we summarize the connection between moments of structure functions 
and parton distributions in different factorization schemes, which in NLO is not unique. For 
the moment, ^2n{Q'^)i of non-singlet part of F2{x,Q'^) we have: 

:F^,!{Q') = c^,!mAr{Q\ (Bl) 

where C^^(Q^) is the Wilson coefficient, and A^^iQ"^) is obtained from the matrix element 
of the local operator taken between nucleon states. Both quantities on the r.h.s. are fac- 
torization scheme dependent. In the MS factorization scheme, for example, the moments of 
parton distributions are related to A^^i^Q"^) simply by: 

In the other widely used factorization scheme, namely the DIS scheme, the Wilson co- 
efficient C2n{Q'^) is absorbed into the definition of the moments of parton distributions 
themselves: 
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(B3) 



Therefore J^^niQ"^) -DIS scheme is defined to have the same form as in the simple 

parton model with free quarks. 

In analogy, the singlet part of the structure function F2{x,Q^) in the DIS scheme must 



also coincide with the simple parton model result |16|: 



(B4) 



Here the Wilson coefficient, as well as the gluonic contribution, are absorbed into the singlet 
quark distribution. However, this procedure does not uniquely fix the gluon distribution. 



For this one enforces momentum conservation, which leads to the relation 16 



^DIS 



MS 



7MS 



(B5) 



To fully define the moments of the gluon distribution one extends this equality to all n. 

Since the anomalous dimensions are in a one-to-one correspondence with the renormal- 
ization constants of local operators, the evolution equations are affected by transformations 
to other factorization schemes. Because the renormalization constants depend on the renor- 
malization scheme chosen to renormalize the local operators (= factorization scheme), the 
anomalous dimensions are factorization scheme dependent. Transforming from the MS to 



the DIS scheme, the anomalous dimensions change according to p3|,H4| 



7 



(1),DIS 



nonsinglet , 



(B6) 



I V 



-'2,n 



singlet. 



From this follow the evolution equations in the DIS scheme. 
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FIG. 1. Ratio of evolution results for qy = uy + dy obtained by calculating moments numer- 
ically and analytically (gf?^). For comparison, two methods for calculating the moments 
numerically (interpolation and histogram) have been applied, using various numbers of points in 
Bjorken-x (given on top of the lines accompanied by a letter indicating the method used). 
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FIG. 2. Comparison of LO and NLO (DIS-scheme) evolution of the total valence distribution 



x{uy + dy)- Input taken from the latest CTEQ fit [11| is evolved from = 100 GeV^ down to 
different values of (given on top of the corresponding lines) . 




35 





36 





37 



